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ABSTRACT 

We present the first cosmological galaxy evolved using the modern smoothed particle 
hydrodynamics (SPH) code GASOLINE2 with superbubble feedback. We show that 
superbubble-driven galactic outflows powered by Type II supernovae alone can produce 
L* galaxies with flat rotation curves with circular velocities ^ 200 km/s, low bulge-to- 
disc ratios, and stellar mass fractions that match observed values from high redshift 
to the present. These features are made possible by the high mass loadings generated 
by the evaporative growth of superbubbles. Outflows are driven extremely effectively 
at high redshift, expelling gas at early times and preventing overproduction of stars 
before z = 2. Centrally concentrated gas in previous simulations has often lead to 
unrealistically high bulge to total ratios and strongly peaked rotation curves. We show 
that supernova-powered superbubbles alone can produce galaxies that agree well with 
observed properties without the need for additional feedback mechanisms or increased 
feedback energy. We present additional results arising from properly modelled hot 
feedback. 

Key words: hydrodynamics - methods: N-body simulation - ISM - galaxies: forma¬ 
tion - galaxies:evolution - cosmology:theory 


1 INTRODUCTION 

The theory of galaxy formation is a cornerstone of modern 
cosmology. ACDM accurately predicts the detailed proper¬ 
ties of the Cosmic Microwave Background, and the forma¬ 
tion of large-scale cosmic structure. Understanding how this 
yields the small-scale properties of the galaxies we see today 
and at high redshift requires a good understanding of the 
physics at play within these individual galaxies: star forma¬ 
tion, feedback, gas accretion, etc. 

The basic theoretical picture of galaxy forma¬ 
tion is now well-established (Rees & Ostriker 1977; 
White & Rees 1978). The complex details are typi¬ 
cally probed through both simulation and semi-analytic 
techniques (Kauffmann et al. 1999; Bower et al. 2006). 
Simulators typically employ large-scale cosmological 
boxes and zoomed in simulations, in which regions sliced 
out of those larger volumes are re-simulated at higher 
resolution (Katz & White 1993; Governato et al. 2004; 
Stinson et al. 2010; Brooks et al. 2011). Recent studies 
(Dekel & Birnboim 2006; Woods et al. 2014) have shown 
that the picture of how gas is fed to galaxies is not simple: 
cold flows can be funneled along filaments in the cosmic 
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web, bypassing the virial shock and directly supplying gas 
to the inner galaxy. Gas accretion and expulsion by stellar 
feedback is fundamental to how galaxies evolve over time, 
determining not only when and how many stars form, but 
the kinematic properties of those stars as well. 

Stellar feedback has a greater impact than merely reg¬ 
ulating star formation by heating the ISM or increasing its 
turbulence. It has long been understood that the cumulative 
effects of multiple supernovae can eject gas from the galaxy, 
powering a galactic outflow or wind (Mathews & Baker 
1971; Larson 1974). Galactic winds are common in high- 
redshift galaxies, and can be seen in the nearby Universe in 
starburst galaxies (See review by Veilleux et al. 2005). The 
detection of both blue-shifted rest-frame UV absorption in 
observations of high redshift galaxies (Weiner et al. 2009; 
Steidel et al. 2010; Martin et al. 2012) along with broad¬ 
ened Ha emission lines (Heckman et al. 1987; Genzel et al. 
2011; Newman et al. 2012) strongly suggests the presence of 
hot, outflowing material surrounding these galaxies. Galac¬ 
tic winds appear to be ubiquitous at high-redshift, and are 
likely a key factor in the evolution of galaxies in the early 
Universe. Muratov et al. (2015) estimated that mass load¬ 
ings for winds peaked at roughly rj = 10 at high 

redshift and that a significant fraction of ejected material 
can escape beyond the virial radius. By ejecting material 
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from a galaxy’s disc, and storing it in the hot, gravitation¬ 
ally bound circumgalactic medium (CGM), star formation 
at high redshift can be strongly regulated, while at the same 
time providing a reservoir of gas that can be accreted at 
late times, ensuring that star formation at low redshift can 
continue (Marasco et ah 2012). How this CGM is created is 
complex, as it is potentially pierced by cold flows of infalling 
material, outflows from the central galaxy, and the accretion 
of dwarf satellites. 

The first generations of cosmological simulations in¬ 
cluding baryonic physics (hydrodynamics, star formation, 
etc.) suffered from a number of serious problems. Nearly 
every simulation, regardless of halo mass, produced far 
too many stars, and tended to form these stars too early 
(Abadi et al. 2003; Governato et al. 2009; Stinson et al. 
2010; Brooks et al. 2011). Simulations of disc galaxies pro¬ 
duced bulge dominated galaxies rather than the thin discs 
we observe. Together with the discrepancies seen between 
large dark-matter only simulations and observations, it ap¬ 
peared that the standard dark energy + cold dark matter 
cosmology ACDM might need to be modified to adequately 
explain the properties of both individuals galaxies and pop¬ 
ulations of galaxies that we see. Scannapieco et al. (2012) 
showed these problems were ubiquitous among codes, lead¬ 
ing simulators to consider a role for stronger feedback. Early 
feedback models (e.g. Katz & White 1993) tended to deposit 
the energy of feedback into a poorly-resolved region of the 
interstellar medium, subjecting it to over-cooling. Stronger 
feedback can be achieved by limiting cooling or increasing 
the total energy. It must also ensure that feedback energy 
couples strongly to the gas. This strong coupling doesn’t 
just heat the ISM, or disrupt only the densest regions, but 
drives outflows that actively remove gas from the disc of 
the galaxy. The Aquila comparison of 13 different simula¬ 
tion codes and subgrid physics models (Scannapieco et al. 
2012) showed that only those cases with strong outflows 
could produce galaxies with realistic star formation histo¬ 
ries. The outflow models used were somewhat ad hoc and 
extremely aggressive. The cases that were capable of eject¬ 
ing sufficient gas from the galaxy disc to moderate the stel¬ 
lar mass of the halo still failed to produce stellar discs and 
largely shut down low redshift star formation. 

Galactic winds have been a part of cosmological galaxy 
simulations for some time (Springel & Hernquist 2003), and 
recent simulations have investigated these winds in detail. 
Angles-Alcazar et al. (2014) showed that these winds can 
dramatically alter the star formation history, kinematics, 
and morphology of galaxies at redshift 2. By explicitly creat¬ 
ing galactic winds with a variety of mass-loadings and wind 
velocities, they showed that strong winds are essential to 
producing the gas-rich, extended, and turbulent discs that 
are typically observed in high redshift star forming galaxies. 
Unfortunately, without simulations to z = 0, interpreting 
exactly how these high-redshift winds impact present-day 
galaxies is difficult. If outflowing material falls back onto 
the galactic disc within a Hubble time, the effects of high- 
redshift winds may be seen in the form of increased star 
formation and inflows at low redshift. 

A key question remains: What processes set 
the mass loading and velocity of these winds? 
Early work tied galactic outflows to the SNe en¬ 
ergy (Springel & Hernquist 2003). Some studies (e.g. 


Murray et al. 2005; Krumholz & Thompson 2013) have 
suggested that radiation pressure is needed to drive suffi¬ 
cient galaxy-scale winds. Others argue that galactic winds 
could be powered by cosmic ray buoyancy (Ipavich 1975; 
Breitschwerdt et al. 1991; Socrates et al. 2008). 

Today, multiple different models of strong feedback have 
begun to produce galaxies with the correct number of stars 
(Aumer et al. 2013), reasonable star formation histories 
(Stinson et al. 2013; Agertz & Kravtsov 2014; Munshi et al. 
2013), and correct morphologies (Guedes et al. 2011; 
Brook et al. 2012; Ghristensen et al. 2014). Unfortunately, 
many of these successes have come at the cost of increasing 
complexity in star formation and feedback methods, crude 
assumptions regarding the physics of the feedback-heated 
gas and somewhat arbitrary increases in the feedback en¬ 
ergy per unit mass in stars. 

In many cases, strong feedback simply means more en¬ 
ergy. Many modern feedback models augment the energy 
input from Type H supernovae (~ 10^^ erg per star above 
8 M©), with that arising from stellar winds, UV ioniza¬ 
tion, supermassive black holes (SMBH), or radiation pres¬ 
sure (Vogelsberger et al. 2013; Agertz & Kravtsov 2014). 
Eeedback models such as Stinson et al. (2013) group these 
sources of energy as early stellar feedback. Since the first 
supernovae occur ~ 4 Myr after star formation, these feed¬ 
back mechanisms begin depositing feedback before SN-only 
methods would. Unfortunately, how much energy from these 
early feedback effects actually couples to the surrounding 
ISM rather than radiating away is highly uncertain. In fact, 
even the coupling of comparatively simpler SN feedback still 
is a matter of some debate. Many new feedback models 
(Agertz et al. 2013; Aumer et al. 2013; Hopkins 2013) treat 
each of these feedback mechanisms explicitly, modelling the 
input of energy and momentum from each component sepa¬ 
rately. 

In addition to increasing the total amount of energy 
deposited by stellar feedback, these models often also in¬ 
clude components designed to prevent energy from feed¬ 
back radiating away (a problem discussed thoroughly in 
Thacker & Gouchman 2000). The energy can be prevented 
from cooling completely for a while (Stinson et al. 2013) or it 
can be initially placed into a non-cooling reservoir that leaks 
back into regular thermal energy (Agertz et al. 2013). Alter¬ 
nately, depositing thermal feedback into a sufficiently small 
mass, allows it to always heat gas to the same high temper¬ 
ature where cooling times are long (Dalla Vecchia & Schaye 
2012). Depositing feedback as kinetic avoids initial radiative 
losses Springel & Hernquist (2003); Agertz et al. (2013); 
Hopkins (2013). 

While these techniques do help solve the overcool¬ 
ing problem, they all come with some drawbacks. Eixed- 
temperature thermal feedback is stochastic, and requires the 
additional free parameter of the feedback temperature. Gool- 
ing shutoffs completely disable radiative losses, where in na¬ 
ture these losses are suppressed in some regions but can re¬ 
main strong in others, depending on the structure of the ISM 
and the clustering of stars. Kinetic feedback is almost always 
paired with a temporary decoupling of hydrodynamic forces 
on feedback-accelerated gas. This is necessary to prevent 
this gas from shock-heating and potentially reintroducing 
the overcooling problem (as shown by Greasey et al. 2011; 
Durier & Dalla Vecchia 2012). Decoupling allows winds to 
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escape. However, this makes it impossible to study the de¬ 
tailed behavior of these winds, and how they interact with 
the ISM. This makes the mass loading an imposed parame¬ 
ter. 

Of these methods, Dalla Vecchia & Schaye (2012) 
showed the interesting result that, for supernovae alone, 
depositing feedback energy into a pre-specified amount of 
mass, without any cooling shutoffs, can give reasonable 
star formation rates and strong galactic outflows in isolated 
galaxies. Keller et al. (2014) argued that what sets this mass 
is thermal conduction in superbubbles, and used that mech¬ 
anism to build a new way of simulating supernovae feedback 
that lacked the resolution dependence, additional complex¬ 
ity, and ad-hoc additions of many current feedback models. 
By focusing on superbubbles, and the evaporation of cold 
gas to determine mass loading, superbubble feedback gives 
realistic gas behavior, and is effective at both regulating star 
formation and driving galactic outflows without introducing 
free parameters. 

The original McMaster Unbiased Galaxy Simulations 
(MUGS; Stinson et al. 2010) showed that with stellar feed¬ 
back, the observed color-magnitude relationship and Tully- 
Fisher relation could be produced in simulated L* galaxies. 
It failed, however, to produce galaxies with the correct stel¬ 
lar mass fraction and star formation history, overproducing 
stars over the entire cosmic history, and grossly overproduc¬ 
ing them at high redshift. It also produced galaxies with 
bulge-to-disc fractions larger than those seen in nature, and 
with sharply peaked rotation curves. Stinson et al. (2013) 
showed that the addition of early stellar feedback could alle¬ 
viate most of these problems. That model does not take into 
account the potentially complex coupling of stellar winds or 
radiation pressure to the surrounding ISM. Instead, a sub¬ 
stantial (fixed) fraction of the stellar bolometric luminosity 
was applied as feedback heating. 

In this paper, we begin first with an overview of the sim¬ 
ulation methods used, both gas microphysics and the star 
formation and feedback techniques. We then present the re¬ 
sults of a suite of 4 simulations using an initial condition 
from the MUGS sample, each generated with different stellar 
feedback models or energetics. The resulting galaxy proper¬ 
ties, as well as the halo evolution are examined over the life¬ 
time of the galaxy. Finally, we discuss how the superbubble- 
driven outflows change with time, and how they ultimately 
result in a realistic galaxy at the present epoch. 


2 METHODS 

These simulations were run using a modern update to 
the SPH code GASOLINE (Wadsley et al. 2004), GASO- 
LINE2. The changes in this new code include a sub-grid 
model for turbulent mixing of metals and energy (Shen et al. 
2010), and a modified pressure force form similar to that 
proposed by Ritchie & Thomas (2001) which is functionally 
equivalent to Hopkins (2013). These changes solve the prob¬ 
lems seen in Kelvin-Helmholtz and blob destruction tests 
with SPH (Agertz et al. 2007). These and other features are 
discussed in Keller et al. (2014). Accurately modelling mix¬ 
ing in multiphase gas is essential for accurately simulating 
the ISM and GGM. 


2.1 Simulations 

Eor this initial study, we have selected the initial conditions 
from one of the original MUGS galaxies. We selected an 
intermediate-mass halo, gl536, allowing us to compare to a 
number of other studies that have examined this particular 
halo (e.g. Stinson et al. 2013; Woods et al. 2014). At z = 0 
this halo has a virial mass ofSxlO^^M© and a spin param¬ 
eter of 0.017. It had its last major merger at z = 2.9. We 
have a gas mass resolution of Mgas = 2.2 x 10^ M©, and use 
a gravitational softening length of e = 312.5 pc. The details 
of how this IG was created can be found in Stinson et al. 
(2010). We choose to focus on a single galaxy for this paper 
as it allows us to make a direct comparison of the impact 
feedback physics makes with relatively little expense. Natu¬ 
rally, this limits our ability to comment on population-wide 
effects. We leave this discussion for a forthcoming paper in 
which we introduce 17 additional L* galaxies. 

We compare four test cases using the same initial 
conditions: no stellar feedback (an absolute lower bound 
for looking at the effects of feedback); superbubble feed¬ 
back (our fiducial case); blastwave feedback (the feedback 
method used in Stinson et al. (2010), first described in 
Stinson et al. 2006); and superbubble feedback with dou¬ 
ble the standard feedback energy (2 x 10^^ erg per SN). 
The high feedback energy case uses more feedback than 
is predicted by Leitherer et al. (1999), but is within the 
range of feedback energies currently being used in cosmo¬ 
logical simulations (Schaye et al. 2015; Agertz & Kravtsov 
2014; Vogelsberger et al. 2013). 

Halos were found in each of the simulations using the 
Amiga Halo Einder (AHE; Knollmann & Knebe 2009). 

2.2 Star Formation and Feedback 

Our suite of simulations use a range of different feedback 
processes. Eor all of the simulations shown in this paper, we 
use a common star formation prescription. Stars are formed 
at a rate proportional to the local free-fall time of gas, such 
that 


Eor each of these simulations, we used an efficiency param¬ 
eter c* = 0.1, the value used by Stinson et al. (2013). Stars 
are allowed to form in a converging flow when gas is cooler 
than 1.5 x 10^ K, and with a density set to that allowed by 
the gravity resolution, p — Mgasj^ — 9.3 cm“^. 

The amount of supernova feedback per unit stellar mass 
is determined using a Ghabrier (2003) IME. With ~ 10^^ erg 
per supernova this gives ~ 10^® ergMQ^. A notable dif¬ 
ference between this simulation and those of Stinson et al. 
(2013) and others is that we do not include early feedback, 
processes such as stellar winds and radiation that can inject 
energy before supernovae occur ~ 4 Myr after the first mas¬ 
sive stars form. A primary role for early stellar feedback is to 
disrupt the densest molecular gas. In these simulations, this 
dense molecular gas cannot be resolved, never being formed 
(and thus it does not form or need to be destroyed). 

The feedback recipe used in our main simulations is 
the superbubble method presented in Keller et al. (2014). 
This method deposits feedback into resolution elements in a 
brief multi-phase state. These particles each have separate 
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specific energies, masses and densities, for the hot bubble 
and surrounding cold ISM, which includes the swept up shell. 
This allows the method to calculate a separate density and 
cooling time with the respective densities for each phase, 
rather than the average density of both phases. Multiphase 
particles are also prevented from forming stars: the average 
temperature of the two phases is essentially always above 
our temperature threshold for star formation. More details 
on this method can be found in Keller et al. (2014). 

The addition of thermal conduction and evaporation in¬ 
troduces some additional time step constraints in order to 
ensure the stability of integration. However, since the ther¬ 
mal conduction rate is capped by the saturation imposed 
by the electron soundspeed, this time step is at worst 1/17 
the Courant time step. In fact, we see an overall speedup 
when using superbubble vs. blastwave feedback, as gas can 
no longer exist in the regime of high density and tempera¬ 
ture (and thus small Courant times). The average compu¬ 
tation time per step was 25% faster using superbubble 
feedback compared to blastwave feedback. These benefits 
are dominant late in the run, as the total amount of dense 
gas becomes larger. We do see some slight increased cost 
early in the run, with a slowdown of ~ 50% before z=4. The 
additional time-step constraints here are similar to the time- 
step constraints required for other methods as well, such as 
decoupled kinetic winds (Springel & Hernquist 2003). 

2.3 Gas Cooling and Physics 

We adopt the same gas cooling physics as a number of 
past simulations using GASOLINE (Stinson et al. 2013; 
Keller et al. 2014). The method for cooling we use here 
was originally presented in Shen et al. (2010). Our sim¬ 
ulations use cooling rates from CLOUDY (Ferland et al. 
2013), and include a redshift-dependent UV background, 
Compton cooling, and primordial and metal cooling from 
10 to 10^ K. This sets these simulations apart from many 
past simulations (Governato et al. 2009; Brook et al. 2011; 
Guedes et al. 2011), which did not include high tempera¬ 
ture metal cooling. We impose an artificial pressure floor, 
using a method described by Robertson & Kravtsov (2008) 
to prevent spurious Jeans fragmentation (as cold, dense gas 
has both Jeans length and Jeans mass below the resolution 
limit of our simulation). We also enforce that the minimum 
SPH smoothing length a particle may have is e/4. This is 
equivalent to a density floor of 400 cm“^ (note that for 
two-phase superbubble particles, this is the maximum mean 
density). These are comparable to the parameter choices in 
another recent re-simulation of the MUGS initial conditions, 
Stinson et al. (2013). 


3 RESULTS 

3.1 Redshift zero Disc Properties 

As can be seen in Figure 1, the galaxy produced with su¬ 
perbubble feedback is disc-dominated and blue. The stel¬ 
lar disc in the case without feedback is nearly non-existent, 
with an ellipsoidal morphology. The superbubble disc ap¬ 
pears somewhat thicker and truncated compared to the disc 
produced with blastwave feedback. The truncation can be 


seen in the reduced stellar scale length of the superbubble 
galaxy, 2.9 kpc vs. 3.8 kpc for the superbubble vs. blast- 
wave galaxies. The thickening appears to be quantitatively 
insignificant, however, as both discs have stellar scale heights 
of 0.9 kpc This may be due to the the disc being disrupted 
at large radii, where the surface densities are lower, allow¬ 
ing superbubbles to grow larger and escape more easily. The 
increased feedback case is interestingly somewhat tilted com¬ 
pared to the rotation of the halo as a whole. 

The gas component of this galaxy can be seen in Fig¬ 
ure 2. The superbubble feedback is clearly much more effec¬ 
tive at ejecting gas from the disc, creating a halo of clumpy 
HI gas. In the doubled feedback case, the disc is heavily dis¬ 
rupted, with no spiral structure visible in either the stellar 
or gas density images. Thus the larger stellar disc is directly 
related to a more extended gaseous disk. 

The rotation curve in Figure 3 shows that this ejection 
also gives a flattened rotation curve without the central peak 
seen in simulations with blastwave or no feedback. This tells 
us that the mass distribution is much less centrally concen¬ 
trated when superbubble feedback is used, more evidence 
that superbubble feedback is effective in preventing bulge 
formation. 

One of the primary problems in older simulations was 
a much larger bulge-to-total fraction based on a kinematic 
decomposition. We adopt the same kinematic decomposi¬ 
tion as the original MUGS study. We calculate, for each star 
within the halo, a circularity parameter, which is simply the 
ratio of the specific angular momentum component perpen¬ 
dicular to the disc (jz) and the specific angular momentum 
for a perfect circular orbit in the same potential that the 
star sees jdrc- As in the original MUGS study, we iden¬ 
tify bulge stars as having jzljdrc < 0.7 and orbital radii 
of < 5 kpc. We find, using these criteria, that our fidu¬ 
cial B/T ratio at z = 0 is 0.09, slightly smaller than the 
Milky Way’s 0.14, and greatly reduced compared to the 
0.55 found in the original MUGS paper, and well within the 
observational constraints from Allen et al. (2006). As these 
numbers, and the distribution of this circularity parameter 
found in Figure 4 show, the vigorous expulsion of central 
gas with superbubble feedback is a powerful bulge preven¬ 
tion and destruction mechanism. Each of the 13 samples 
from the Aquila project (Scannapieco et al. 2012) suffered 
from either from a bulge-dominated stellar circularity pro¬ 
file, or from a massively peaked (or in some cases simply too 
high) rotation curve. Even the three cases that managed 
to produce a reasonable stellar mass fraction (each of these 
generated strong outflows, either through SMBH feedback 
or wind decoupling) still failed to produce disc-dominated 
galaxies. Of their sample, only four cases had more than 
40% of the stellar mass in disc stars, but all of these cases 
massively overproduced stars. 

The ejection of gas is evident if we look at the baryon 
content of the interior part of the halo, where the disc re¬ 
sides, in this case, simply material within O.lRvir,- With 
superbubble feedback, the baryon fraction of this inner re¬ 
gion is 0.30, reduced from 0.37 without feedback, and 0.35 
with blastwave feedback. This means that nearly 20% of the 
baryons that would be available to form stars have been 
blown out of the galaxy disc, roughly twice the amount that 
was removed by the older feedback model. The baryon deficit 
in the superbubble is comparable to the total stellar mass 
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No Feedback Blastwave Superbubble Superbubble, 2X FB 



Figure 1. Mock stellar image of each galaxy at 2 ; = 0. The top row shows the galaxies edge-on, while the bottom shows them face-on. 
The no-feedback case shows very little disc, while the blastwave feedback case has a thin disc with a prominent bulge. The superbubble 
galaxy appears nearly bulgeless, composed only of a truncated disc. Note also that the superbubble galaxy is bluer than the blastwave 
or no-feedback galaxies, evidence of a younger stellar population. 



Figure 2. HI column density for the four test cases at z=0. The no feedback case has exhausted nearly all of the gas within the disc, 
leaving only diffuse wisps. The two superbubble cases (especially with doubled supernova energy) have lofted a large amount of gas out 
of the disc, and some entrained dense clumps can be seen in this outflowing material. The superbubble gas discs are much more flocculent 
than the blastwave disc, and where the feedback energy is doubled, the disc is strongly disrupted. 


Simulation 


^vir -^gas -^star -^hary ,0.1vir ^bulge ^disc 


No Feedback 

7.94 X 10“ 

6.21 X IQio 

7.42 X lOio 

Blastwave Feedback 

7.99 X 10“ 

8.05 X 10i° 

5.11 X 10i° 

Superbubble Feedback 

8.02 X 10“ 

1.19 X 10“ 

1.84 X lOio 

Superbubble, 2 x 10^^erg/SN 

7.96 X 10“ 

1.20 X 10“ 

1.09 X 10i° 


8.09 X IQio 
7.94 X 
6.25 X IQio 
4.92 X lO^o 


9.34 X 10^ 
8.09 X 10® 
1.14 X 10^ 
5.96 X 10® 


2.88 X lO^o 
1.75 X lO^o 
1.15 X lO^o 
6.98 X 10® 


Table 1. Halo components at z = 0. MfUsc cind are the determined by the angular momentum of the stars (as detailed below). 

All masses are in Mq. 
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Figure 3. Rotation curves for each of our test cases. As is clear, 
the superbubble cases have much lower central concentrations, 
giving a a rotation curve that rises to a flat 200 km/s. The peaked 
rotation curves in the cases with blastwave or no feedback are 
a result of their failure to remove low angular momentum gas 
at high redshift, giving bulge-dominated, centrally concentrated 
galaxies. 



Jz/Jcu 


Figure 4. Histogram of stellar orbit circularity. The bimodal 
distribution shows the bulge component, with a circularity near 
0, and the disc component, with a circularity near 1. Note that 
without feedback, this bimodality disappears, as the galaxy mor¬ 
phology becomes spheroidal. The values here are normalized so 
that each curve has a total integral of unity. As is clear, very few 
stars produced in a galaxy with superbubble feedback have low 
angular momentum, giving us a disc-dominated galaxy. 

of the galaxy (the fiducial case has a 1.5 x 10^° M© deficit 
in baryons (compared to the no-feedback case), and a total 
stellar mass of 1.8 x 10^°M©). 

Of what remain, only 29% of the baryonic mass in the 
superbubble galaxy disc is in stars, compared to 63% in the 
blastwave galaxy, and 89% in the no feedback galaxy. The 
basic properties of each halo can be found in table 1. 



Figure 5. Stellar mass growth as a function of the total halo 
mass. The grey band and black dashed curve show mean val¬ 
ues and uncertainties from the abundance matching results 
of Behroozi et al. (2013). The points show values at redshifts 
3,2,1,0.5, and 0.1. The stars show the final values at redshift 0. For 
essentially the entire evolution of the halo, both the no feedback 
galaxy and the blastwave galaxy overproduce stars. 


54 3 2 1 0.5 0.1 0 



Time {Gyr) 

Figure 6. Stellar mass fraction as a function of time/redshift. As 
in figure 5, only the flducial superbubble run produces stellar mass 
fractions within observed uncertainties for the entire history of the 
halo. The grey band is as in figure 5, and here we also show, in 
the hatched region with the dotted curve, values and uncertainties 
from Moster et al. (2013). The dip that can be observed in the no 
feedback case at high redshift come from the growth of the halo, 
as it accretes dark matter dominated dwarves. 


3.2 Halo Evolution and Star Formation 

Ejecting baryons from the disc is essential to producing 
galaxies that fit the M^/Mhaio relationship predicted by 
Behroozi et al. (2013), Moster et al. (2013) and others. As 
Figures 5 and 6 show, the galaxies either lacking feedback 
or using blastwave feedback alone fail to regulate star for¬ 
mation, diverging from the expected abundance matching 
relation at z ~ 3 for the blastwave, and before z = 5 with¬ 
out feedback, giving galaxies that lie above the abundance- 
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z 


54 3 2 1 0.5 0.1 0 



Time (Gyr) 


Figure 7. Star formation rates for all 4 simulations. The low 
star formation rate seen after 2 ; = 1 in the no feedback case is 
due simply to the high z star formation consuming most of the 
available gas. 


5 4 3 2 



Figure 8. Star formation rate up to z = 2 for our superbubble 
galaxy. As is clear, when sampled on short (7 Myr) timescales, the 
star formation rate in the superbubble galaxy is characterized by 
strong bursts, despite its quiescent behavior averaged over longer 
timescales. 


matching relationships over nearly their entire history and 
mass evolution. With superbubble feedback, this galaxy lies 
within the range of observed stellar masses over its entire 
evolution, although on the low side of this range at low red- 
shift. Arbitrarily increasing feedback energies, despite hav¬ 
ing reasonable star formation rates near z = 0 , under pro¬ 
duces stars for most cosmic history, giving stellar masses 
lower than predicted by abundance matching. Stellar and 
total masses were calculated for the region inside Rvir- As 
the halo grows over time, brief reductions in the stellar mass 
fraction can be seen for the no-feedback and blastwave feed¬ 
back halos in Figure 6 , simply by the accretion of gas and 
dark matter that has not yet made it to the galaxy disc. 
Once these mergers complete, the stellar fraction rate leaps 
up once again because of the new supply of gas. The relative 
flatness in this figure at low redshift is reflected in the rel¬ 
ative flatness in the star formation rate, shown in figure 7. 
This is simply a side effect of the massive star formation 
rates at high redshift: the bulk of star forming gas has been 
consumed, and thus star formation has slowed. 

The star formation rate, in bins of 150 Myr, is shown 
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Figure 9. Galactic wind mass-loading and outflow rates for each 
feedback model. As is clear, superbubble feedback powers winds 
with significantly higher mass-loadings, and, despite forming less 
than half as many stars the blastwave galaxy, results in more 
total outflowing gas. This is key to the regulatory power of super¬ 
bubbles. Despite the peak in outflow rate in the blastwave case 
between 2 ; = 4 — 1.5, the massively increased star formation rate 
seen in the the smoothed SFR during this time (bottom plot) 
means that the mass-loadings never go above 5, and baryon ex¬ 
pulsion is inefficient. 


in Figure 7. With no feedback, or blastwave feedback, star 
formation is extremely vigorous before z — 2^ slowing at 
low redshift as gas available for star formation is exhausted. 
With superbubble feedback, star formation is relatively con¬ 
stant over nearly the entire history of the halo, with only a 
gradual increase towards z = 0. This apparent quiescence 
is a function of the large time window over which we are 
averaging our star formation rate. In Figure 8 , you can see 
clearly that despite having an average star formation rate of 
~ 1 Moyr“^, this is punctuated by bursts of star formation 
of as much as ~ 10 M 0 yr“^. Very similar results were seen 
in Muratov et al. (2015), where bursts of star formation are 
followed by peaks in the outflow rates. This burstiness is a 
important if stellar feedback is to drive galactic winds: to 
generate large superbubbles, star formation must be clus¬ 
tered in both space and time. Whether these bursts of star 
formation are able to effectively remove baryons from the 
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Figure 10. Net baryonic accretion (stars and gas) for each test 
case. It is clear that superbubble feedback expels baryons much 
more effectively than the blastwave feedback model (although 
even weak feedback does help somewhat compared to no feed¬ 
back whatsoever). 





Figure 11. Inflow and outflow rates for superbubble (red) and 
blastwave (blue) run. The solid lines show the total inflow rate, 
while the dashed and dotted lines show rates for cold and hot gas 
(above/below 10^ K) respectively. 


disc ultimately depends on the mass loading of the winds 
that they drive. 


3.3 Outflow Analysis 

Star formation in the disc drives hot, fast-moving outflows 
from the central star forming regions. These outflows have 
temperatures of~ 2 xlO^Kas they leave the disc, and also 
entrain some cold material with them. Typically outflow ve¬ 
locities are a few hundred km/s, less than the escape velocity 
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Figure 12. In all four test cases, including the fiducial su¬ 
perbubble run shown above by the solid curve, and the no¬ 
feedback case shown by the dotted curve, the vast majority of 
low angular momentum material, with specific angular momen¬ 
tum jz < 500 kpc km/s, is accreted by 2; = 2. The accretion of 
high angular momentum material, with jz > 1500 kpc km/s rises 
to peak at z = 2, and continues to accrete ~ 2 M0/yr at 2; = 0. 
As can be seen here, the net flux of low angular momentum ma¬ 
terial at high redshift is suppressed, while the accretion of high 
angular momentum material at low redshift is actually enhanced 
by superbubble feedback. 
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Figure 13. High-redshift winds preferentially remove low angular 
momentum material. As can be seen here, significantly more low 
angular momentum material is removed by superbubble feedback 
than by the weaker blastwave feedback, but the amount of ejected 
high-angular momentum material is roughly equal. 


of the halo, but sufficient to propel gas to large radii before 
it begins to fall in again. 

We calculated inflow and outflow rates and velocities 
by examining particles within a spherical shell, with inner 
radius ^.IRvir and outer radius Rvir (giving a shell of thick¬ 
ness 0.9Rvir)- We use the p — 299pcrit definition for Rvir- 
Particles with Vr <9 are inflowing, while those with Vr > 9 
are outflowing. Within this shell, the mass flux is determined 
simply as: 


M= ^ 

Gshell 


MjVj 

0.9Rvir 


( 2 ) 


As figure 9 shows, the outflow behavior of this galaxy is 
fundamentally diherent when superbubble feedback is used. 
The wind mass loading factor 77 = Moutfiow /a measure 
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Figure 14. Phase diagram of gas in the halo at ^ = 0. The 
hot, coronal gas with temperatures above 10^ K contains roughly 
40% of the gas mass in the halo. This gas is a mixture of virial- 
shocked pristine gas, which has never been within the interior 
O.lRyir, and outflowing material ejected from the galaxy. This 
buoyant, high-entropy gas exits the galaxy at high temperature, 
and cools adiabatically as it rises through the CGM. The three 
colored contours show the major components of this hot coronal 
gas. Pristine gas, never falling within the interior, is shown in blue. 
Young superbubbles, not yet having broken out of the galaxy, are 
shown in red. Outflowing material that was once inside the galaxy 
but is now cooling adiabatically as it rises through the CGM can 
be seen in orange. The contours show the 1.5%, 1%, and 0.5% 
levels for the total mass in each cut. 


of how efficiently feedback is generating outflows, is roughly 
an order of magnitude larger during the bursts of star for¬ 
mation from z = 4 — 2, expelling a large amount of gas from 
the disc early on. The resulting suppression of high-redshift 
star formation is key to obtaining the correct stellar mass 
relation as the halo grows, is unmistakable in figure 5 and 6 . 
Evidence of these outflows can be seen in the z = 0 column 
density images in figure 2. High-latitude neutral hydrogen 
from both ejected and infalling material can be seen in the 
two superbubble cases, but is totally absent otherwise. The 
smaller amount seen around the high feedback energy case 
is simply due to this gas being expelled further. The rea¬ 
son superbubbles are so effective at removing mass from the 
galaxy and thus regulating the stellar mass depends both 
on the higher outflow rates as well as the lower star forma¬ 
tion rates needed to drive these outflows, leading to much 
larger mass loading factors than the blastwave model can 
achieve. This is clearly shown in figure 9, as the top panel 
is essentially the middle panel divided by the bottom panel. 
The outflow rates driven by superbubbles is at most ~ 2 
times greater than the blastwave model, but since this is 
driven by a comparatively tiny amount of star formation, 
the massloading, 77 , is nearly 10 times greater at high red- 
shift. The factor of ~ 2 higher outflow rate is what we would 


expect from the roughly ~ 2 larger baryon depletion from 
the superbubble vs. blastwave cases. 

The effectiveness of the superbubble-driven galactic 
winds can be clearly seen in figure 10 , where the net ac¬ 
cretion of stars and gas is greatly reduced by the use of 
superbubble feedback, resulting at z = 1 in a net reduc¬ 
tion in the total baryonic mass inside O.bRvir of ~ 14% vs. 
blastwave feedback, and ~ 30% vs. no feedback whatsoever 
(along with the removal of baryons from the disc discussed 
earlier. As the outflow rate drops towards z = 0 , much of 
this ejected material falls back onto the disc, and the out¬ 
flows transition to a fountain mode, actually increasing the 
inflow rate seen compared to the blastwave feedback case 
in figure 11. This increase in the accretion of gas fuels an 
equivalent increase in star formation, as is seen in figures 7 
and 9. 

Figure 12 shows that the total accreted gas switches 
from being dominated by low angular momentum mate¬ 
rial to being primarily high angular momentum material 
at z ~ 2. Because rj is high (~ 10) during this period, 
low-angular momentum material is preferentially ejected 
from the galaxy without forming a significant number of 
stars. These results agree quite well with the results of 
Muratov et al. (2015), which also found L* progenitors at 
z > 2 had wind mass loadings 77 ~ 10. The preferential ejec¬ 
tion of low angular momentum gas is clear in figure 13. The 
small amount of high angular momentum material accreted 
before z = 1 is ejected at an essentially equal rate with or 
without superbubble-driven winds, but low angular momen¬ 
tum material is propelled into winds at roughly twice the 
rate when superbubble feedback is taken into account. 

The phase diagram in figure 14 shows that, as expected, 
the galaxy contains both a hot halo (in which nearly half of 
the gas mass can be found), and contains no gas in regions 
of short cooling time. This behavior was shown previously in 
the isolated galaxy simulations of Keller et al. (2014). The 
hot halo gas can be divided into three major components. 
The coolest component is pristine, virial shocked material 
that has not ever been accreted (never passed within ^.IRvir 
of the halo center). The hottest gas is actually not yet in 
the corona, but is the interior of young superbubbles still 
within the galaxy disc. As this material leaves the disc, it 
cools adiabatically in the lower pressure environment of the 
halo, as can be seen in from the adiabatic path it takes 
through the phase diagram. Detailed analysis of this halo 
gas, especially in comparison to observations like those of 
Steidel et al. (2010), will be explored in a future study of a 
larger sample of L* galaxies. 


4 DISCUSSION 

Past work, especially the Aquila comparison 
(Scannapieco et al. 2012 ), has shown that feedback mecha¬ 
nisms which do not remove a large fraction of baryons from 
galaxy discs fail to produce realistic spiral galaxies. Galaxies 
simulated without such feedback show more spheroidal 
stellar distributions, older stellar populations, and more 
total stellar mass than is observed in nature. Those cases 
that do manage to expel enough gas to produce reasonable 
stellar mass fractions still fail to produce the correct 
stellar kinematics, failing to produce the small bulges 
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characteristic of so many spiral galaxies. These models also 
rely on mechanisms other than stellar feedback (SMBH 
heating) or on major numerical contrivances (hydrodynamic 
decoupling, etc.). 

Superbubbles offer a way out of this bind. As was shown 
in (Keller et al. 2014), evaporation in superbubbles natu¬ 
rally mass load winds with temperatures of 10^ AT, a 
process that is set by self-regulating thermal conduction. 
This means that an optimal amount of material is heated 
above the virial temperature (1.3 x 10® K for a 8 x 10^^ Mq 
halo). This feedback-heated and highly buoyant hot gas mi¬ 
grates out of the disc, cooling adiabatically while it rises, 
as is clearly seen in figure 14. In fact, as figure 6 shows, 
there may even be room to reduce the feedback energy below 
the fiducial 10®^ erg/SNe, while still producing reasonable 
stellar mass fractions. On top of yielding physically-realistic 
galaxies, superbubble feedback is in fact slightly less compu¬ 
tationally expensive than blastwave feedback for these sim¬ 
ulations, due in part to its elimination of unphysical high 
density, high temperature gas. 

4.1 High-Redshift Outflows Determine Galaxy 
Properties 

A major failure of other feedback models is the production 
of too many stars, too early. Superbubble feedback prevents 
this by efficiently removing gas from the star forming disc at 
high redshifts. The high mass loadings from z = 2 —4 means 
more mass is being expelled by significantly fewer stars. 
Mass loadings much larger than unity have been observed in 
both Lyman Break Galaxies at high redshift (Pettini et al. 
2002) and in local dwarf starburts (Martin et al. 2002). 

Even if a galaxy has the correct stellar mass fraction 
at z = 0, forming too much of your stellar mass at high 
redshift is problematic for a number of reasons. First, it is 
known from abundance matching as well as the observations 
of stellar ages that most stars in L* galaxies are formed fairly 
late (e.g. van Dokkum et al. (2013) found ~ 90% of stars in 
Milky Way like galaxies formed at z < 2.5). Secondly, form¬ 
ing these stars early on means they will be formed primarily 
in smaller halos that are subsequently accreted, as well as in 
main halo at small radii (as low angular momentum material 
is accreted first, as seen in figure 12). This results in a stel¬ 
lar distribution too spheroidal and centrally concentrated, 
as was seen in many older simulations, and here in the ro¬ 
tation curves in figure 3 for the blastwave and no feedback 
test cases. 

Because superbubbles drive strong outflows at high red¬ 
shift, the galaxy is able to preferentially remove low angu¬ 
lar momentum gas, as indicated by figure 13. This same 
mechanism was seen in simulations of dwarf galaxies by 
Brook et al. (2011, 2012). We would suggest that this is a 
universal requirement for producing galaxies with low cen¬ 
tral concentrations and small bulges. As Binney et al. (2001) 
argued, if this gas were to remain within the disc, it would be 
ultimately form the massive bulge that is seen in our blast- 
wave and no feedback cases, and in numerous other past sim¬ 
ulations. If only a third of the baryons removed in the disc of 
our fiducial case were to instead form bulge stars, it would 
raise the B/T ratio of the disc to 0.3, making our disc much 
more bulge dominated than the Milky Way, and putting 
our results in conflict with Allen et al. (2006). This mech¬ 


anism (strong, high-redshift outflows), and the fact that it 
produces strong suppression of star formation at high red¬ 
shift as well as a small stellar bulge, agrees well with the 
results of Sales et al. (2012). That study of ~ 100 GIMIC 
(Grain et al. 2009) galaxies found that those with the small¬ 
est B/T ratio also had the largest fraction of stars formed 
late in the galaxy’s history. Strong outflows acting to reg¬ 
ulate high redshift star-format ion while preserving fuel for 
late star formation naturally leads to this situation. 

The effectiveness of superbubble feedback may help to 
prevent the growth of a massive stellar bulge via a second 
mechanism as well. Fall & Efstathiou (1980) showed that 
the dissipational collapse of hot gas within an extended 
dark halo can produce galaxies with thin discs. However, 
Gole et al. (2000) showed that the collisionless processes in¬ 
volved in galaxy mergers can cancel out the net rotation of 
a galaxy disc and lead to a spheroidal stellar distribution. 
With superbubble feedback, small dwarves convert very lit¬ 
tle of their baryonic mass into stars, allowing them to con¬ 
tribute their gas to the primary halo, which then collapses 
dissipationally. 

The heavily mass loaded winds driven by superbubble 
feedback are not only key to suppressing high-z star forma¬ 
tion, but also to providing enough gas at low redshift to con¬ 
tinue star formation. The outflows produced by superbubble 
feedback can easily escape from the lower-mass progenitors 
of L* galaxies. Thus it is able to naturally transition from 
the violent, wind-driving high redshift mode to a more qui¬ 
escent phase as the galactic gravitational potential and gas 
surface densities in the disc increase (see the increase in in¬ 
flow in figure 11. These factors suppress outflows, and allow 
star formation to increase slowly towards z = 0. 

As the disc is assembled, the disc surface density above 
the superbubble increases. Thus the hot gas must push past 
more material and will also entrain more material along with 
it as shown by the clumps seen in figure 2). This slows the 
outflows compared to those seen at high redshift and mov¬ 
ing the galaxy towards a fountain mode. The gas is kicked to 
relatively low altitudes and then rains back down onto the 
galactic disc. This effect is probably sensitive to resolution. 
Our small-scale experiments show that venting of superbub¬ 
bles is enhanced by a more porous ISM (see Keller et al. 
2014 and also Nath & Shchekinov 2013). Poor resolution 
suppresses this porosity and increases numerical dissipation. 
Thus better resolution might reduce the puffy appearance in 
the HI column density and allow strong galactic outflows to¬ 
wards z = 0. 

4.2 Additional Feedback Mechanisms 

We have shown that, even at moderate resolution, thermal 
supernova feedback is sufficient to build a realistic L* galaxy, 
provided a complete physical feedback model like that of 
Keller et al. (2014) is used. Thus this work firmly establishes 
what supernovae and superbubbles can do. 

It is still likely that for higher-mass halos (> 10^^ M©), 
supernovae alone will not be suffficent to regulate star for¬ 
mation and produce the drop in star formation efficiency 
seen in Moster et al. (2013); Behroozi et al. (2013). We do 
not include feedback from SMBH, which is likely to be im¬ 
portant for higher-mass galaxies. 

Our resolution limits the formation of dense structures 
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in the ISM. The lack of resolved dense gas means that feed¬ 
back processes involved in disrupting molecular clouds (UV 
photodissociation, radiation pressure, etc.) is largely irrele¬ 
vant for these simulations. In fact, since much of the energy 
from early stellar feedback is consumed in the disruption of 
the densest clouds in simulations such as this, the addition 
of this energy may be unrealistic, effectively double counting 
energy that would have been used to disrupt clouds. Further¬ 
more, high-resolution studies of individual molecular clouds 
have found that ionizing radiation, despite disrupting these 
clouds, ultimately only imparts < 0.1% of the total radia¬ 
tive luminosity to the gas in the cloud and the surrounding 
medium in the form of additional thermal and kinetic energy 
(Dale et al. 2005; Gendelev & Krumholz 2012; Walch et al. 
2012)). Thus, applying even as little as 1% of the radiative 
luminosity as a source of feedback in a simulation that does 
not resolve structure within molecular clouds is likely mas¬ 
sively overestimating the impact on disc scales. 

In sufficiently high-resolution simulations it may be 
necessary to include small scale feedback mechanisms in 
order to disrupt clouds before the first SN, allowing 
the to explode in a lower density environment since, as 
Mac Low & McCray (1988) showed, the cooling time for 
superbubbles scales sub-linearly with the inverse density, 
tR oc . Rogers & Pittard (2013) showed that the dis¬ 

ruption of high-density clouds by stellar winds prior to the 
first supernovae can allow as much as 99% of the hot SN 
ejecta to escape into the surrounding warm ISM, helping 
to promote the growth of superbubbles that then vent from 
the galactic disc. In the current work, the densities near the 
superbubbles are modest and this preprocessing is not re¬ 
quired. 


5 CONCLUSION 

We have shown that supernova feedback alone, with a com¬ 
plete physical superbubble model, is capable of producing 
an L* galaxy that falls within observational constraints. Su¬ 
perbubbles are a physical mechanism for producing ab initio 
galactic winds, that ultimately allow for the formation of a 
galaxy with a realistic star formation history and a negligible 
stellar bulge. 

The key results with respect to galaxy formation are as 
follows: 

• Strong outffows at high redshift are essential to regu¬ 
lating star formation over the total halo history. The vast 
majority of stars in L* galaxies form after z ^ 2. Unless gas 
is removed from galaxies at high redshift by feedback pro¬ 
cesses, it will rapidly form stars, yielding discs that are too 
massive and red at z = 0. 

• Outflows are important for producing the correct disc 
kinematics and preventing the formation of excessive bulges. 
Low angular momentum gas is accreted first as galaxies 
form, and the pooling of this gas at the center of galaxies 
can lead to galaxies with sharply peaked rotation curves and 
unrealistically large bulges, often containing the majority of 
stars within the galaxy. 

• Galaxies simulated without feedback, or with disabled 
cooling feedback models fail to expel enough of this gas early 
on, and result in bulge-dominated galaxies with unrealistic 
stellar mass fractions. The superbubble model, on the other 
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hand, produces strong outflows that ultimately yield realistic 
galaxies. 

• Superbubble feedback naturally yields the sort of out¬ 
ffows that are required for L* galaxies. The mass-loading 
& velocity of the winds are set by the hydrodynamics and 
evaporative mixing, unlike other methods where these values 
are free parameters. 

• Superbubble feedback produces high-redshift outflows 
that preferentially remove low angular momentum gas, and 
in so doing, prevents the formation of massive bulges and 
the associated strongly peaked rotation curves. It does this 
without also expelling additional high-angular momentum 
gas, allowing the stellar disc to form while arresting the for¬ 
mation of a bulge. 

• These advantages come without any significant addi¬ 
tional computational expense, and may in fact be less costly 
than alternative feedback models. 

Superbubbles are effective up to at least L*. Beyond 
L*, we anticipate important roles for SMBH feedback or 
potentially radiation pressure (see e.g. Hopkins 2013). In 
the subsequent MUGS2 series of papers, we will show how 
these effects extend to larger and smaller halos forming in a 
range of environments. 
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